
library(dplyr)
library(tidyr)
library(readr)
library(stringr)
library(fixest)

test_data_path = 'cfr test'
mdl = 'NOmatching' # change to 'matching' for matching model

D = vector(mode = 'list',length = 2)
names(D) = c('math','reading')
for (sbj in names(D)){
  D[[sbj]] = file.path(test_data_path,paste0('tv_',sbj,'_',mdl,'.csv')) %>%
    read_csv() %>%
    arrange(districtid,schoolid,schgrade,year) %>%
    group_by(districtid,schoolid,schgrade) %>%
    mutate(lag_tv_2yr_f = lag(tv_2yr_f),
           lag_score = lag(score)) %>%
    mutate(d_score = score - lag_score,
           d_va = tv_2yr_l - lag_tv_2yr_f) %>%
    select(districtid,schoolid,schgrade,year,n,d_score,d_va)
}
D = bind_rows(D, .id = 'sbj') %>%
  mutate(cohort = year - schgrade)

m1 = feols(d_score ~ d_va | schgrade, data = D, weight = ~n, cluster = ~districtid^schoolid^cohort)
nobs(m1)
round(coeftable(m1),3)
round(confint(m1),3)

m2 = feols(d_score ~ d_va | schgrade + schoolid^year, data = D, weight = ~n, cluster = ~districtid^schoolid^cohort)
nobs(m2)
round(coeftable(m2),3)
round(confint(m2),3)

m3 = feols(d_score ~ d_va | schgrade + schoolid^year^sbj, data = D, weight = ~n, cluster = ~districtid^schoolid^cohort)
nobs(m3)
round(coeftable(m3),3)
round(confint(m3),3)







